A first attempt at implementing automatic MC statistical uncertainties, à la Combine.

I followed the Combine implementation outlined here.

Combine handles automatic MC stats by using text2workspace.py, which parses the datacard and finds all TH1 objects and creates a RooRealVar object for all param objects in the card. These objects get attached a Gaussian constraint and are treated as a nuisance parameter. Notably, this only works for RooWorkspace objects that do not contain RooDataHists. This disqualifies 2DAlphabet workspaces from using automcstats, since we use 2D RooDataHists.

I think that enabling something like what Combine does using params in the datacard would necessitate a very difficult re-write of the Combine code, so this implementation uses a more naive approach. During workspace creation, when building the OrganizedHist object, 2DAlphabet will now loop over all nominal background histograms in each region (signal can be included as well, though I don't see why one would want to). The algorithm follows Combine's implementation:

1. Sum the yields $n_i$ and uncertainties $e_i$ of each background process $i$ in bin (bx, by) of the 2D histogram.
2. Calculate $n_\text{tot} = \sum_{i\in \text{bkg}} n_i$, and $e_\text{tot} = \sqrt{\sum_{i\in \text{bkg}} e_i^2}$
3. If $e_\text{tot}= 0 $, the bin is skipped and no Up/Down $\pm 1\sigma$ TH2s are created
4. Calculate the effective number of unweighted events $n_\text{tot}^\text{eff} = n_\text{tot}^2 / e_\text{tot}^2$
5. If $n_\text{tot}^\text{eff} \leq n^\text{threshold}$, separate uncertainties are created per-process
   • NOTE: Combine creates a Poisson-constrained parameter. We can't really do this with 2DA, so instead we just create $\pm 1\sigma$ templates (where $-1\sigma$ is bounded below at 0), and then Combine automatically assigns these as nuisance parameters with Gaussian constraints.
6. If $n_\text{tot}^\text{eff} > n^\text{threshold}$, then we create a single Gaussian-constrained Barlow-Beeston-lite template that will scale the total yield in the bin.
   • NOTE: I do this by ensuring that every process in this bin shares the same nuisance parameter (histogram) name, so that every background process is controlled by the same parameter in Combine. See below for some reasoning.

The resulting histograms are identical to the nominal background histograms except for the given bin, which is shifted up or down by some scaling factor governed by the number of events.

Case 1: $n_\text{tot}^\text{eff} > n^\text{threshold}$

In this case, we are creating a single nuisance parameter that scales the total yield in the bin. To accomplish this, I use the following logic:

Each process $i$ in this bin gets its nominal yield $n_i$ scaled by $1\pm \text{r}$, where the "down" variation is bounded at 0. For this case, $r = e_\text{tot} / n_\text{tot}$. Thus, in this bin, process $i$ contributes $n_i * (1 + \nu * r)$, where $\nu$ is the value of the (Gaussian-constrained) nuisance parameter. For all processes in this bin, we then obtain:

$$ \begin{align} n_i^\pm &= \sum_{i\in \text{bkg}} n_i * (1 \pm \nu * r) \\ &= (1 \pm \nu * r) * \sum_{i\in \text{bkg}} n_i \\ &= (1 \pm \nu * e_\text{tot} / n_\text{tot}) * n_\text{tot} \\ n_i^\pm &= n_\text{tot} \pm \nu * e_\text{tot} \end{align} $$

Since $e_\text{tot} = \sqrt{\sum_i e_i^2}$, we get the combined MC stat uncertainty in this given bin. And, since $\nu$ is shared for all processes in this bin, a single nuisance parameter controls the yields of all processes in this bin.

Usage

The user does not have to do anything special, I've enabled 2DA to use "autoMCstats" by default. The method can be turned off, its threshold modified, or signals included in the calculation, by modifying the JSON:

    "OPTIONS": {
        "mcstats": true,
        "mcstats_threshold": 10, 
        "mcstats_include_signal": false
    }

Other than that, no change needs to be made to the user's 2DAlphabet python script or JSON config file. Note that this will produce at least N * R * 2 new histograms, where N is the number of processes, R is the number of regions (pass, fail, etc), and 2 is for up/down. There may be more histograms if the BB-lite threshold is not met in a given bin.

Testing

I've tested this with some simple dummy histograms simulating $X \to H(b\bar{b}) Y(b\bar{b})$ in a Fail and Pass region and 2D histograms of data, a single signal, and Zbb and Hbb backgrounds. The testing workspace can be recreated with the following files:

test.tar.gz

The tarball contains:

• an input ROOT file containing the data, Zbb+Hbb backgrounds, and one signal. The backgrounds and signal are subject to two shape uncertainties.
• a JSON file for 2DAlphabet defining the pass and fail regions, all processes, and the two shape systematics. The JSON also contains an OPTIONS section that customizes the mcstats options
• a python script to create the 2DAlphabet workspace with a few different transfer function parameterizations

Testing results

By default, I have 2DAlphabet give verbose output when creating the mcstat nuisances. The output is more or less copied entirely from combine. Here's the example output when creating the workspace:

(twoD-env) [ammitra@cmslpc331 dummy_test]$ python create_workspace.py 
Doing GLOBAL variable replacement in input json...
Running with options: Namespace(verbosity=0, overwrite=False, debugDraw=False, blindedPlots=[], blindedFit=[], mcstats=True, mcstats_threshold=10, mcstats_include_signal=False, haddSignals=True, plotTitles=False, plotTemplateComparisons=False, plotPrefitSigInFitB=False, plotEvtsPerUnit=False, year=1)

================================================================
Analyzing bin errors for region: fail
Poisson cut-off: 10
Processes summed: Zbb Hbb
Processes excluded from sums: XHY-1800-100
================================================================
Bin                 Contents           Error  Notes
----------------------------------------------------------------
fail (1, 1)         6.000000        2.449490  total sum
fail (1, 1)         6.000000        2.449490  excluding marked processes
fail (1, 1)         6.000000        2.449490  Effective events, alpha=1.000000
  => N_eff <= 10 : per-process (Poisson approximated as gaussian)
  ----------------------------------------------------------
    Zbb                      6.000000       2.449490
                             6.000000       2.449490  Effective events, alpha=1.000000
      => 'mcstat_fail_Zbb_bx1_by1' [1.00,0.00,3.86] to be gaussian constrained (width=0.408248)
  ----------------------------------------------------------
    Hbb                      0.000000       0.000000
      => Content or error is zero, ignore
----------------------------------------------------------------
fail (1, 2)        51.000000        7.141428  total sum
fail (1, 2)        51.000000        7.141428  excluding marked processes
fail (1, 2)        51.000000        7.141428  Effective events, alpha=1.000000
  => N_eff > 10 : shared gaussian shape 'mcstat_fail_bx1_by2' (rel=0.140028) shared by [Zbb, Hbb]
----------------------------------------------------------------
...

One can already see two examples of the MC stat functionality working for cases where (1) the $n_\text{tot}^\text{eff}$ is above threshold and (2) it is below the threshold. In the former case a shared parameter is created and in the latter separate parameters are created per-process.

One can also see that signal is excluded from this calculation and from MC statistical uncertainties, as per the JSON options.

We can also look at the post-fit results for the fits with and without MC statistical uncertainties, using a constant transfer function for the QCD estimate. All other fit options are the same between the two cases. In the former, the pulls are noticeably smaller, indicating that the fit is making use of the additional freedom provided by the MC statistical uncertainty parameters.

mX projection:

mY projection: